Basics of BigDFT: N2 molecule as example

This is a simple notebook that shows how to execute a simple calculation with BigDFT. You will learn how to manipulate basic DFT objects from a python script.

How to perform a first run with default parameters

For this tutorial, we will run just a simple N2 molecule, taken from the BigDFT database (available here).

from BigDFT.Systems import System
from BigDFT.Fragments import Fragment
from BigDFT.IO import XYZReader

N2 = System()
with XYZReader("N2") as ifile:
    N2["N:0"] = Fragment([next(ifile)])
    N2["N:1"] = Fragment([next(ifile)])

Every system object has a posinp representation, which is a dictionary representation of the geometry. We will use this representation later.

print(N2.get_posinp())

{'positions': [{'N': [0.0, 0.0, 0.5488], 'frag': ['N', '0']}, {'N': [0.0, 0.0, -0.5488], 'frag': ['N', '1']}], 'units': 'angstroem', 'cell': [inf, inf, inf]}

We now run the DFT calculation (with the default input parameters) of this molecule. This can be done by instantiating a class from the Calculators module. Then we call the run method of that calculator with the posinp of N2 as a parameter.

from BigDFT import Calculators as C
from os import environ
environ['I_MPI_FABRICS']='shm'
study = C.SystemCalculator(verbose=False) #Create a calculator
log = study.run(posinp=N2.get_posinp(),name="N2") #run the code

The run method of the BigDFT.Calculators.SystemCalculator class shows in the standard output which is the command that is invoked. Then, an instance of the Logfile class is returned. Such class can be used (we will see after) to extract the information about the electronic structure of the system.

Now we can retrieve important informations on the run. See the examples below:

log.energy #this value is in Ha

-19.884615273244822

log.evals[0].tolist() # the eigenvalues in Ha ([0] stands for the first K-point, here meaningless)

[[-1.041353673208,
  -0.4926440819324,
  -0.4357815638016,
  -0.4357812404171,
  -0.3818323179544]]

log.get_dos().plot() #the density of states

<matplotlib.axes._subplots.AxesSubplot at 0x153189d2b040>

The tutorial DoS Manipulation can be followed to learn how to plot Density of States from a Logfile

How to use the SystemCalculator instance

Before managing different calculations at once, let us look closer at the SystemCalculator that was presented above.

study = C.SystemCalculator(verbose=False,omp=1,mpi_run='mpirun -np 2')

As you can see, initiating an instance allows to set the computational parameters such as the OpenMP and MPI parallelisation. In this case, a single thread is used while two processes are running in parallel, therefore reducing the computation time. The global options of the runner (or calculator) can then be accessed by

study.global_options()

{'omp': '1',
 'mpi_run': 'mpirun -np 2',
 'dry_run': False,
 'skip': False,
 'verbose': False}

Global options can be also added and removed as follows

study.update_global_options(new_option = 'value')
study.global_options()

{'omp': '1',
 'mpi_run': 'mpirun -np 2',
 'dry_run': False,
 'skip': False,
 'verbose': False,
 'new_option': 'value'}

Similarly, global options can also be removed

study.pop_global_option('new_option')
study.global_options()

{'omp': '1',
 'mpi_run': 'mpirun -np 2',
 'dry_run': False,
 'skip': False,
 'verbose': False}

How to modify the input parameters from the default ones

To specify non-default input parameters to the code, we should employ a Inputfile object. One can choose a naming prefix for a run, which enables to classify the runs which are performed, and ease the subsequent postprocessing.

Methods of the class can be employed for modifying the input parameters. For instance, the XC functional can be chosen via the set_xc method. Imagine for example that you are interested in visualizing the wavefunctions output of the calculation. To do that, the suitable method of the class instance should be called. Create a new calculation set by using, for instance, the LDA prefix.

from BigDFT import Inputfiles as I

inp = I.Inputfile()
inp.set_xc('LDA')
inp.write_orbitals_on_disk()

The wavefunctions will be present at the end of calculation, by indicating the value of the key orbitals in the output dictionary.

#Run the code with the LDA prefix
LDA = study.run(name="LDA",input=inp,posinp=N2.get_posinp())

When using a naming scheme, the output files are placed in a directory called data-LDA. In our LDA example, the wavefunctions of the system can thus be found in the data-LDA directory:

!ls data-LDA

time-LDA.yaml                     wavefunction-k001-NR.bin.b000004
wavefunction-k001-NR.bin.b000001  wavefunction-k001-NR.bin.b000005
wavefunction-k001-NR.bin.b000002  wavefunction-rhoij.bin
wavefunction-k001-NR.bin.b000003

Here 001 means the first K-point (meaningless in this case), N stands for non spin-polarized, R for real part and the remaining number is the orbital ID. Post-processing of these files is done in the another tutorial.

In the same spirit, another calculation can be done with different parameters. Imagine we want to perform a Hartree-Fock calculation. We should simply change the XC functional used. However, we also have to specify the pseudopotential.

inp.set_xc('HF')
inp['psppar.N']={'Pseudopotential XC': 1} #this set the default PSP for LDA
HF = study.run(name="HF",input=inp) #Run the code with the name scheme HF

Without specifying the PSP, an error will occur. The same error is specified at the end of the log file, and also in debug/bigdft-err-0.yaml:

 Additional Info:
   The pseudopotential parameter file "psppar.N" is lacking, and no registered pseudo found
   for "N"

This is because the pseudopotential is assigned by default in the code only for LDA and PBE XC approximations. Another alternative is to specify the internal PSP that might be used, taking from the default database of BigDFT. This can be done as follows:

inp['psppar.N']={'Pseudopotential XC': 1} #here 1 stands for LDA as per the XC codes

When possible, care should be taken in choosing a pseudopotential which has been generated with the same XC approximation used. Unfortunately, at present HGH data are only available for semilocal functionals. For example, the same exercise as follows could have been done with Hybrid XC functionals, like for example PBE0 (ixc=-406). In the case of Hartree-Fock calculations, using semilocal functionals generally yield accurate results (see Physical Review B 37.5 (1988): 2674).

In BigDFT, XC functionals are specified using the built in named functionals, or using the LibXC codes.

inp.set_xc('PBE0')
PBE0 = study.run(name="PBE0",input=inp) #Run the code with the name scheme PBE0

The variables LDA, HF, and PBE0 contains all information about the calculation. This is a class Logfile which simplify considerably the extraction of parameters of the associated output file log-LDA.yaml. If we simply type:

print(LDA)

- Atom types:
  - N
- cell: Free BC
- number_of_orbitals: 5
- posinp_file: LDA.yaml
- XC_parameter: -20
- grid_spacing:
  - 0.45
  - 0.45
  - 0.45
- spin_polarization: 1
- total_magn_moment: 0
- system_charge: 0
- rmult:
  - 5.0
  - 8.0
- dipole:
  - -0.0005134858
  - -0.0005134858
  - -0.000559708
- energy: -19.88461527324464
- fermi_level: -0.3818323179544
- forcemax: 0.01156401391398
- forcemax_cv: 0.0
- gnrm_cv: 0.0001
- nat: 2
- symmetry: disabled
- No. of KS orbitals:
  - 5

We display some information about the LDA calculation. For instance, we can extract the eigenvalues of the Hamiltonian i.e. the DOS (Density of States):

LDA.evals[0].tolist()

[[-1.041353673208,
  -0.4926440819324,
  -0.4357815638016,
  -0.4357812404171,
  -0.3818323179544]]

Note that LDA.log is the python dictionary associated to the full output file. This tutorial explains the basic methods.

print(PBE0)

- Atom types:
  - N
- cell: Free BC
- number_of_orbitals: 5
- posinp_file: PBE0.yaml
- XC_parameter: -406
- grid_spacing:
  - 0.45
  - 0.45
  - 0.45
- spin_polarization: 1
- total_magn_moment: 0
- system_charge: 0
- rmult:
  - 5.0
  - 8.0
- dipole:
  - -0.0004219248
  - -0.0004219248
  - -0.0004512606
- energy: -19.971244525934665
- fermi_level: -0.4512290316272
- forcemax: 0.05514199331393
- forcemax_cv: 0.0
- gnrm_cv: 0.0001
- nat: 2
- symmetry: disabled
- No. of KS orbitals:
  - 5

The following exercise (and its solution) gives some clues about it.

The overall structure of files in the disk

In the disk, after that the run is performed, different files are created:

• input_minimal.yaml which contains the minimal set of input variables to perform the run;

• log.yaml which contains the log of the run with all calculated quantities;

• time.yaml and forces_posinp.yaml which we will see later.

For its I/O, BigDFT uses the yaml format. If you look at the input_minimal.yaml file, you can see:

  #---------------------------------------------------------------------- Minimal input file
  #This file indicates the minimal set of input variables which has to be given to perform
  #the run. The code would produce the same output if this file is used as input.
 posinp:
   units: angstroem
   positions:
   - N: [0.0, 0.0, 0.0]
   - N: [0.0, 0.0, 1.114989995956421]
   properties:
     format: xyz
     source: posinp.xyz

In this case, only the atomic positions are indicated using a yaml format.

However, there is a one-to-one correspondence between a yaml file and a python dictionary. For this reason we will create the input parameters of this runs from dictionaries. See the useful yaml online parser page to understand the correspondence.

In the log file log.yaml, BigDFT automatically displays all the input parameters used for the calculations.

The parameters which were not explicitly given are set to a default value, except the atomic positions of course, which have to be given. As we did not specified input files here, our run is a single-point LDA calculation, without spin-polarisation.

Exercise

Compare the values of the HOMO and HOMO-1 eigenvalues for the LDA and the HF run. Change the values of the hgrid and crmult to find the converged values. Note that a thorough description of those two (essential) parameters is provided later on when studying a CH4 molecule.

Note that, both in the LDA and in the HF calculation, a norm-conserving PSP is used.

The results can be compared to all-electron calculations, done with different basis sets, from references (units are eV)

1. S. Hamel et al. J. Electron Spectrospcopy and Related Phenomena 123 (2002) 345-363 and (2) P. Politzer, F. Abu-Awwad, Theor. Chem. Acc. (1998), 99, 83-87:

eigenvalues	LDA(1)	HF(1)	HF(2)	(Exp.)
3σg	10.36	17.25	17.31	(15.60)
1πu	11.84	16.71	17.02	(16.98)
2σu	13.41	21.25	21.08	(18.78)

The results depends, of course, on the precision chosen for the calculation, and of the presence of the pseudopotential. As it is well-known, the pseudopotential appoximation is however much less severe than the approximation induced by typical XC functionals. We might see that, even in the HF case, the presence of a LDA-based pseudopotential (of rather good quality) does not alter so much the results.

Here you can find the values from BigDFT calculation using a very good precision (hgrid=0.3, crmult=7.0). Note that 1 Ha=27.21138386 eV.

eigenvalues	LDA	HF
3σg	10.40	17.32
1πu	11.75	16.62
2σu	13.52	21.30

How much do these values differ from the calculation with default parameters? Do they converge to a given value? What is the correlation for the N2 molecule in (PSP) LDA?

See here the solution.